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Abstract. At nonzero chemical potential the numerical sign problem in lattice field theory limits the use 
of standard algorithms based on importance sampling. Complex Langevin dynamics provides a possible 
solution, but it has to be applied with care. In this review, we first summarise our current understanding 
of the approach, combining analytical and numerical insight. In the second part we study SL(7V, C) gauge 
cooling, which was introduced recently as a tool to control complex Langevin dynamics in nonabelian gauge 
theories. We present new results in Polyakov chain models and in QCD with heavy quarks and compare 
various adaptive cooling implementations. 

PACS. 12.38.Gc Lattice QCD calculations 
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1 Introduction 

. One of the outstanding questions in the theory of the 
■ strong interactions. Quantum Chromodynamics, is the de- 
, termination of its phase diagram in the plane of tem- 
' perature and baryon chemical potential. Nonperturbative 
I studies using lattice QCD are hindered by the presence of 
t the complex fermion determinant and the resulting sign 
problem: since the Boltzmann weight in the partition func- 
tion is complex, standard numerical approaches relying 
I on importance sampling are only of limited use [T] . As a 
I result various alternative approaches in numerical lattice 
field theory have been developed during the past years; a 
comprehensive overview can be found in Ref. (2, . 

In this review we focus on the complex Langevin (CL) 
' method, in which a complexified configuration space is ex- 
plored stochastically and the presence of a positive weight 
is not required [31I3]. In fact, due to the complex nature 
of the Boltzmann weight, the dynamics naturally takes 
place in this larger complexified configuration space, and 
it is precisely this extension that allows for a solution of 
the sign problem. While it has been successfully demon- 
strated that the sign problem can be solved in a number of 
theories, including those in which the sign problem is se- 
vere [5,ifiji7,J,i success is not guaranteed [9, 10, ll,12, 13 . 14t 
[T51I16) . The first aim of this contribution is to summarise 
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why problems may arise and how they can be detected in 
practice. Recently, definite progress has been made in non- 
abelian gauge theories, by employing gauge cooling during 
the complex Langevin evolution |17| . The second aim of 
this paper is to discuss this in some detail and present new 
results comparing various gauge cooling implementations. 

This contribution is organised as follows: in the fol- 
lowing section, we contrast complex and real Langevin dy- 
namics, and remind the reader why the standard proof for 
the justification of real Langevin dynamics is no longer ap- 
plicable. A justification following a different route is given 
in Sec. [31 we identify criteria for correctness and empha- 
sise the importance of understanding the distribution in 
the complexified space. This is illustrated in Sec. HI where 
two three-dimensional models are contrasted: the SU(3) 
spin model and the XY model, both at nonzero chemi- 
cal potential. While for the former CL works in the entire 
phase diagram, for the latter this is only true in the or- 
dered part of the phase diagram. In Sec.[5]we come to non- 
abelian dynamics and explain how gauge cooling can be 
used to control the distribution in the complexified space, 
using Polyakov chain models and QCD with heavy quarks 
for illustration. It is shown that a judicious choice of cool- 
ing implementations leads to well-controlled dynamics in 
SL(iV, C). Ref. [2] reviews CL with an emphasis on aspects 
complementary to the ones discussed here. 



2 Real versus complex Langevin dynamics 

In order to indicate the difference in theoretical status be- 
tween real and complex Langevin dynamics, we start with 
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a brief reminder, emphasising why the "standard proofs" 
of correctness are no longer available in the case of a com- 
plex action. The results in this section are well known [T5j . 

For simplicity, we consider one degree of freedom x, 
but the expressions can easily be generalized. We start 
with the partition function 

Z^jdxe-^^''\ S'(a;)eR. (2.1) 

The corresponding Langevin equation reads 

X = -d.,S{x) + T^, (2.2) 

where the random noise rj is Gaussian distributed, with 
{r]{t)ri{t')) = 2S{t—t'). The evolution of expectation values 
of observables 0{x) in Langevin time is represented, after 
averaging over the noise, by the evolution of the associated 
distribution p{x,t), according to 

(O)p(t) = J dxp(x,t)0{x). (2.3) 

Given the Langevin equation for x{t), it follows that p{x, t) 
satisfies a Fokker-Planck (FP) equation, 

p{x,t)=d^{d^ + S'{x))p{x,t). (2.4) 

This FP equation has a stationary solution p{x) ^ e~^^'^\ 
which is the first prerequisite. Moreover, it turns out that 
this solution is typically reached exponentially fast. In or- 
der to see this, write p{x,t) = ^{x,t)e~^^''^^; it then fol- 
lows that ^{x,t) satisfies a Schrodinger-like equation, 

^ix,t) = -HFPij{x,t), (2.5) 

with the FP Hamiltonian, 

Hfp - Q«Q, (2.6) 

in which 

= -d^ + \s'{x), Q = d., + ^S'{x). (2.7) 

Since — Q'' , the FP Hamiltonian is semi-positive def- 
inite: if goes to infinity as |a;| does, it has a discrete 
spectrum and a nondegenerate ground state (cf. Refs. [121 
[T5]). Let us denote the eigenvalues of iJpp by A > 0, with 
A = corresponding to 

Q^{x)^0 ^ V(a;) ~ e^s^C^), (2.8) 

then 

^{x, t) = coe-5^^(") + ^ cxe-^* ^ coe-^^("\ (2.9) 

A>0 

and the correct distribution is reached exponentially fast. 

In the case of a complex action, S{x) E C, and assum- 
ing that S can be analytically continued into the complex 



plane, one can attempt to follow the same line of reason- 
ing. Replacing x by z to emphasise that the dynamics is 
now complex, the Langevin equation 

i = -9^5(2) +7? (2.10) 

and FP equation 

p{z,t)^dAd, + S'{z))p{z,t) (2.11) 

can still be written down. However, since in this case 
p{z, t) is complex-valued, it is not a probability distribu- 
tion. Moreover, the associated FP Hamiltonian (j2.6p is no 
longer semi-positive definite, since ^ . 

In order to discuss the (real and positive) probability 
distribution in the complex plane, we write explicitly x — > 
z — X + iy and consider the complex Langevin equations, 

x^K^ + r]R, y^Ky + rji, (2.12) 

with the drift terms 

= -Red,S{z), Ky = -lmd,S{z). (2.13) 

Here we momentarily introduced "complex" noise, writing 
■q = r]n + im ^MM- Starting from {r]{t)r]{t')) = 2S{t-t'), 
it follows that the real and imaginary components of the 
noise satisfy 

{m{t)vi{t')) =mS{t-t'), (2.14) 

imitW)) =0, 

with ATj^ - A^i = 1 and A'^i > 0. The evolution of the 
expectation value of holomorphic observables 0{x) is now 
represented by 

(O)p(t) - J dxdy P{x, y- t)0{x + ly), (2.15) 

where the associated distribution P{x, y; t) satisfies the 
(real) FP equation 

P{x,y;t) = [d^{Nnd,-K,) 

+dy{Nidy-Ky)]P{x,y-t). (2.16) 

Although this equation is superficially of the same form as 
Eq. (j2.4p . with simply twice as many degrees of freedom, 
it is in fact much harder to solve. In particular there are 
no generic solutions known, even in the stationary case. 
Some recent solutions in special cases can be found in 
Refs. POIITSII^ . Moreover, most likely it is not possible 
to express the dynamics in terms of a semi-positive defi- 
nite FP Hamiltonian, so that convergence to a stationary 
solution is not guaranteed. The best one may hope is a 
Hamiltonian whose spectrum is in the right half plane, 
but this should be sufficient. Nevertheless, in many cases 
of interest a stationary distribution emerges when the CL 
equations are solved numerically: one goal of this review 
is to indicate how properties of this distribution relate to 
the reliability and justification of the approach. 
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We end this section with a description of how the CL 
equations appear in lattice field theory. Consider first a 
real scalar field (j)ax, where a denotes additional inter- 
nal degrees of freedom. After discretising Langevin time 
with time step e and using the lowest-order discretisation 
scheme, the Langevin equation takes the form 



xin + l) = 



with the drift 



(t)ax{n) + eKax{n) + V^Vaxin), 



Kax = -SS[<l>]/S<l>a 



(2.17) 



(2.18) 



and the noise (we only consider real noise from now on, 
Ni = 0) satisfying {riax{n)r]a'x'{n')) = 2(5aa'4x"5„„'. When 
the drift is complex, the field is complexified, (pax — 



hlx- 



iipaxi above. 



For gauge or matrix theories, with gauge group SU(A^), 
the dynamics is written for group elements, living on links, 
as [31111] 



Ux,,,{n + 1) = Rx,u{n) Ux,y{ri), 



with 



Rx,. 



exp 



axu I 



(2.19) 



(2.20) 
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Here Xa (a = l,...iV^ — 1) are the traceless, hermitian 
Gell-Mann matrices, normalized as TrAaA^, ~ '^Sab, and 
the drift is determined by 



(2.21) 



where the action S[U] may contain a contribution from the 
fermion determinant as Sp = — IndetM. The differential 
operator D is defined as 



a=0 



(2.22) 



In the case that the action is complex, it is easy to see 
that the drift if 7^ if t, and hence i?ti? ^ 1, ahhough 
its determinant deti? = 1. In SU(iV) theories, one finds 
therefore that U G SL(iV,C) during the CL process. 

Finally, we note that in practice it may be necessary 
to choose the stepsize e adaptively, to avoid instabilities 
due to runaway solutions [22j . In general, discretisations of 
the Langevin equation result in finite stepsize corrections, 
which are linear in e for the Euler scheme given above. This 
can be improved by employing higher-order discretisation 
schemes ,23„24.8,. 



3 Distributions in the complex plane 

A crucial role is played by the distribution P{x,y) in the 
enlarged configuration space (we drop the t dependence to 
indicate that from now on we only consider the stationary 
distribution). While it is not guaranteed that a station- 
ary solution exists in general, in many interesting cases 
we find that it does and that it can be constructed by 



Fig. 1. U(l) one-link model: scatter plots in the complex plane, 
for ^ = 0.1 and 1, at fixed /? = 1 and k = 1/2 [5]. 



brute force, namely by solving the CL equations and col- 
lecting the coordinates visited during the stochastic pro- 
cess. Properties of the distribution can then be studied 
using histograms and scatter plots. This method is well 
suited for zero-dimensional models; when more degrees of 
freedom are involved, it is useful to integrate out most 
coordinates and focus on the ones of interest. 

To illustrate this, we consider a simple U(l) one-link 
model ,5 , with the partition function 



Z= f da;e'''=°"^detM(a;; 

J —7T 



m), 



(3.1) 



with 



det M{x; /i) = 1 + k cos(a:: — i/i). 



(3.2) 

Its structure is inspired by QCD in the hopping expansion, 
with a complex "determinant" , satisfying [det M(x; /i)]* = 
det M{x, —fj,*). When /Lt = 0, the determinant is real and 
the Langevin process takes place on the real axis, y = 
(we take |k| < 1, so that the determinant is positive when 
fj, = 0). Increasing fi results in the trajectories presented 
in Figs. [U [5] Shown here are CL trajectories, indicated 
by the small (blue) dots in the x — y plane (note that the 
dots are separated by 500 updates). The arrows indicate 
the direction of the classical forces, while the larger (black) 
dots are classical fixed points, where the drift terms van- 
ish, Kx = Ky = 0. A small increase in fi moves the dis- 
tribution slightly into the complex plane, see Fig. [T] (top). 
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Fig. 2. As in Fig. [J for ^ = 2 and 5. 

Increasing [i, even further, as in Fig. [2j moves the distri- 
bution further out, but it remains locahsed and bound to 
the attractive fixed point at x = 0. Note that the distribu- 
tion is periodic in the x direction and highly asymmetric 
in the y direction. Expectation values of observables ob- 
tained with CL agree with the exact results, which can be 
computed in this case [S]. 

However, in most cases of interest exact results are not 
available and CL has to be justified without recourse to a 
comparison. For this a good understanding of the proper- 
ties of the distribution is necessary. We found that it is es- 
sential for the mathematical justification of the approach 
that the distribution is well localised in the complexified 
configuration space and in particular shows a fast decay 
in the imaginary directions. The argument goes along the 
following lines [T51[T5] . We consider the two expectation 
values, 



pit) 



dz p{z, t)0{z), 



Pit] 



(3.3) 

+ (3.4) 



in which p{z,t) and P{x,y;t) were introduced above and 
satisfy their respective FP equations, (12. lip and (12. 16^ . 
The goal is to show that 



(3.5) 



for all times or possibly only in the limit that t — ^ oo, start- 
ing from identical initial conditions (although the initial 



condition dependence is expected to drop out). This can 
indeed be shown by moving the time evolution from the 
densities to the observables, making use of the Cauchy- 
Riemann equations and performing integration by parts, 
assuming that no boundary terms at infinity are produced. 
Crucial is therefore the decay at ?/ — >■ ±oo of the product 
P{x, y)0{x + iy) for a suitable class of observables. In 
Ref. [T^ we showed concretely for a toy model how the 
failure of the CL method is associated with insufficient 
decay of P{x, y)0{x + iy). 

In Ref. [15j we also identified certain consistency con- 
ditions which are necessary for the CL method to give 
correct results. They involve the complex "Langevin op- 
erator" 



L = [d,-{d,s{z))]d,, 



(3.6) 



which governs the evolution of holomorphic observables; 
the consistency conditions are 



(LO) = 



(3.7) 



for all observables O. Here the expectation value is taken 
with respect to the equilibrium distribution P{x,y), or 
equivalently after averaging over the noise. Formally these 
conditions just express the stationarity of the equilibrium 
expectation values, but actually their derivation involves 
integrations by part and therefore requires sufficient de- 
cay of the equilibrium distribution. If the decay is insuffi- 
cient, the conditions may fail either by being nonzero, by 
diverging, or by being ill-defined (reflected in huge fluc- 
tuations). While these conditions are only necessary, in 
principle they become sufficient if they hold for a suffi- 
ciently large set (a basis in some sense) of observables 
[T5] . They were also shown to be formally equivalent to 
the Schwinger-Dyson equations. In Ref. [TS] the strength 
of these conditions, applied to a few simple observables, 
was tested and generally we found that they are successful 
in weeding out incorrect results. 

Another intuitive insight relates to the stability of the 
real manifold (i.e. j/ = 0) under small complex fluctua- 
tions (25; . In the example discussed above, the stationary 
distribution at /i = is given by P(x, y; p = 0) = p{x; ^ = 
0)S{y). In Fig. [1] (top) it can be seen that turning on /z, 
and thereby introducing a nonzero force in the imaginary 
direction, does not drastically change the distribution: it 
only slightly widens in the y direction. The distribution 
at /i ^ is therefore connected to the one at /i = 0. Al- 
though this condition is neither sufficient nor necessary 
(there could be another, much wider distribution which 
gives the correct result), it helps in interpreting results in 
practice, as we will demonstrate below. 

Finally, we note that when complex noise is used (iVj > 
0), there is additional diffusion in the imaginary direction. 
This has the effect of widening the distribution. We found 
that in most cases this has a detrimental effect on the re- 
sults of CL [T^ITS] . In general complex noise is therefore 
not recommended, although in some cases a small amount 
of complex noise {Ni <C 1) can be used to regulate numer- 
ical solutions of the FP equation (I2.16P [13] . 
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Fig. 3. SU(3) spin model: (Tr(;7 + !7~^)/2) as a function of 
/^^, for several values of /3 on a 10^ lattice j8j. 

4 SU(3) spin model vs XY model 

In this section we apply these ideas to contrast two three- 
dimensional models, the SU(3) spin model and the XY 
model, both at nonzero chemical potential. It turns out 
that CL works very well in the SU(3) spin model [S], but 
that it fails in part of the phase diagram of the XY model 
[14]. Understanding the difference yields considerable in- 
sight into the approach [23] . We note that both can also be 
solved with worldline/flux methods, which is sign-problem 
free, so that "exact" results are available [2611271128] . 

The three-dimensional SU(3) spin model is an effec- 
tive Polyakov loop model for QCD with heavy quarks, 
obtained using a combined strong-coupling/hopping ex- 
pansion. The action reads 

S^^^Y. C^=.Tr V'l^ + Tr [/"^Tr V^+o\ 

-h KTr [/, + e-^Tr U~^] , (4.1) 

X 

where the JT^'s are SU(3) matrices, living on a three- 
dimensional lattice with fi = sites, and the addi- 
tional sum in the first line is over the nearest neighbours. 
Static quarks are represented by Polyakov loops and cou- 
ple to the chemical potential in the usual way, resulting 
in a complex action, S*{^) = S{—ii*). This model has 
been solved with CL in the classic papers [2^150] : here 
we show results from Ref. [5]. The flux representation 
is discussed in Refs. [271128) , while a recent mean- field 
study appears in Ref. [21] ■ The model is part of a fam- 
ily of strong-coupling/heavy-quark effective models, con- 
structed in Refs. [521133] . We note that in the latter CL 
has been used to study the phase structure. 

The action for the three-dimensional XY model is 
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Fig. 4. XY model: action density {S) /O as a. function of ^'^ in 
the ordered phase at /3 = 0.7 (top) and the disordered phase at 
j3 = 0.3 (bottom) on an 8^ lattice. The dotted line at ^ = ivr/S 
indicates the Roberge- Weiss transition at imaginary ^ |14) . 



Here the degrees of freedom are U(l) phases, Ux = 6**^=". 
The chemical potential couples to the conserved Noether 
charge and the action is complex, as above. This model has 
been studied using a worldline formulation in Ref. :26 and 
more recently in Ref. [33] . Here we discuss the results from 
Ref. T^. Both models have a similar phase structure: a 
disordered phase at small /3 and and an ordered phase at 
large /3 and/or fx. In the SU(3) spin model, the transition 
weakens as /.t is increased and turns into a crossover. The 
critical couplings at ^ = are /3c = 0.133 — 0.137 (depend- 
ing on h) in the SU(3) spin model [5 ] ! ^ and /3c = 0.45421 
in the XY model 35 . 

Out of the number of criteria we have developed to 
justify CL, we consider here (the lack of) analyticity in 

around /i^ = and changes in the distribution in the 
complexified configuration space as is increased. Note 
that negative jj^ corresponds to imaginary chemical po- 
tential, for which real Langevin dynamics can be used. 
Fig. 13 shows (Tr([/ + C/"^)/2) as a function of in the 
SU(3) spin model, for several values of (3. It is clear that 
the transition from /i^ < (obtained with real Langevin) 
to > (obtained with complex Langevin) is smooth 
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Fig. 5. XY model: normalised difference between the CL and 
the worldline results for the action density, AS, in the (3 — 
fj, plane. The phase boundary between the ordered and the 
disordered phase is indicated with the black line [14| . 
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Fig. 6. XY model: width of the distribution in the imaginary 
direction, {(zi<jf>i)^), as a function of ^ for several j3 values ,14 . 

and that the observable is analytic in /i^, as it should be. 
This interpretation has been confirnied with simulations 
using the flux representation [28) . 

This should be contrasted with what happens in the 
XY model. In Fig. 2] we show the expectation value of 
the action density around fi^ ^ 0, for j3 = 0.7 in the or- 
dered phase (top) and (3 — 0.3 in the disordered phase 
(bottom). While in the ordered phase the observable is 
smooth in /i^, in the disordered phase we observe a clear 
nonanalyticity, which is interpreted as a failure of CL. This 
conclusion is further supported by a comparison with re- 
sults obtained in the world-line formulation at positive : 
again we see agreement in the ordered phase but disagree- 
ment in the disordered phase. We also note that the data 
point at /i^ = in the lower line of the bottom graph 
has been obtained using complex initial conditions for the 
Langevin process. The resulting disagreement with the ex- 
pected result indicates that in this case the real manifold 
is unstable and that the distribution at /i ^ is not con- 
nected to the one at = (see below). 

The failure in the XY model turns out to be closely 
related to its phase structure, as demonstrated in Fig. [SJ 



What is shown here is the normalised difference between 
the CL result and the worldline result for the action den- 
sity, 

^S^{SU^A^_ (4.3) 

The phase boundary, indicated with the black line, is taken 
from Ref. [13]. We conclude that incorrect results are ob- 
tained in the disordered/transition region, while CL is 
seen to work in the ordered phase. As an aside, we re- 
mark that these findings are independent of the volume 
used and therefore independent of the severity of the sign 
problem. 

The question is how the mismatch seen in Figs. U (bot- 
tom) and [5] manifests itself in the distribution in the com- 
plexified space, -P[0R, 0i] . Since this is a high-dimensional 
example, it is not possible to plot the distribution directly, 
as in Figs.[Tl[2j However, we can determine how the width 
of the distribution in the imaginary direction, 

{{Acj^if) = {cpD - (4.4) 

depends on /i and /3. This is shown in Fig. [S] At /i = 0, 
((Zi0i)^) — provided that real initial conditions for the 
Langevin process are used. At large /3 the width increases 
smoothly from as /i is increased. At smaller /? on the 
other hand, the width jumps discontinuously. Note that 
we used complex initial conditions, introducing complex- 
ity even at /i = 0. Alternatively, one can obtain this dis- 
continuity using an infinitesimal but nonzero /x. We also 
note that at large /i, deep in the ordered phase, the width 
decreases and the dynamics becomes effectively real, sig- 
nalling that the sign problem is effectively absent [25, . We 
conclude therefore that in the disordered phase the distri- 
bution P[(/'R, at /i ~ is not connected to the distri- 
bution yo[0] at ^ = 0. The distribution P[(j)B., </>!] is broad 
and not sufficiently localised, leading to a breakdown of 
the validity of CL dynamics. 

The difference between the SU(3) spin model and the 
XY model can be understood from an analysis of clas- 
sical flow diagrams, as in Figs. [U [2l after reducing them 
to effective one-link models [25j . The crucial difference be- 
tween the (abelian) XY model and the (nonabelian) SU(3) 
model is the presence of the reduced Haar measure in the 
latter. It is exactly this contribution that stabilises the dy- 
namics, since it is purely restoring and directed towards 
attractive fixed points on the real manifold. On the other 
hand, in the XY model there is no restoring component 
and the real manifold is unstable against small complex 
fluctuations |25) . 

5 Gauge cooling 

The lessons from the previous sections are that the dis- 
tribution in the complexified space should be sufficiently 
localised and that it is desirable that turning on the chem- 
ical potential leads to a smooth deformation of the initial 
distribution into the complexified configuration space. We 
now apply these insights to SU(Ar) gauge theories. As de- 
scribed at the end of Sec. [21 during a CL process links 
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U originally in SU(7V) will naturally take values in the 
noncompact group SL{N,C), in which [/^ ^ U~'^. Be- 
fore making the analytic continuation, we therefore first 
express the action and observables in terms of the holo- 
morphic variables U and U~^, rather than in terms of . 
To determine to which extent SL(A'', C) is explored and 
the distribution is localised in this larger group, a mea- 
sure needs to be defined that signifies the distance from 
the real submanifold, i.e. from SU(iV). There are a number 
of possibilities, here we consider the "distance" [5] 

cf=lTr(C/C/t-]l)>0, (5.1) 
and the "unitarity norm" 

norm = Tr {UU^ - l)^ > 0. (5.2) 

In both cases equality is reached for unitary matrices. The 
inequality in Eq. (|5.ip follows from the polar decomposi- 
tion oi U e SL{N,C), which states that U = PV, with 
V e SU(A^) and P = P'l' a positive-definite matrix with 
determinant det P = 1 . 

We demonstrate the approach using Polyakov chain 
models: Ni links Uk form a one-dimensional chain, with 
periodic boundary conditions, and the partition function 
reads 

Z= / l[dU,e-''^^l (5.3) 

fc=i 

The action, with complex coupling /?, reads 

S=-^TT{UiU2...UNi) (5.4) 

for SU(2) and 

S = -/3iTr(C/i . . . Un,) - p2MUj^^ . . . U^^) (5.5) 
for SU(3), with 

/3i(m) = /3 + Ke'^, 

/32(/i) = r +«;e-'^ = [/3i(-M*)]*- (5.6) 

The structure of these models is motivated by the effec- 
tive one-link models encountered in QCD in the heavy- 
quark/strong-coupling limit [2S]. In the case of SU(2), 
there is no sign problem at nonzero chemical potential, 
but a complex-action problem appears at nonzero theta 
angle or in real time. The coupling in Eq. (15. 4p is chosen 
to be complex to niimick this. 

As in four-dimensional SU(iV) gauge theory there is a 
notion of gauge symmetry: a link starting at site k trans- 
forms as 

Uk ^ nkUk^-l^, = e^^""", (5.7) 

where here and below the sum is over a = 1, ■ . • , N"^ ~ 1. 
Of course, these models are gauge-equivalent to one-link 
models (i.e. with A^^ = 1) and exact results can be eas- 
ily computed by integrating over the reduced Haar mea- 
sure [2S]. However, increasing offers the possibility to 



study how the performance of the CL algorithm depends 
on many gauge degrees of freedom, mimicking the four- 
dimensional case. 

During the CL simulation, a distribution P[[/, C/^] is ef- 
fectively being sampled (we explicitly indicated here that 
P depends both on U and U^). In the case that the out- 
come of the CL simulation is incorrect, we may now exploit 
the gauge freedom to attempt to change the distribution 
P[U, C/t] in SL(A^, C). Note that an SL(A^, C) gauge trans- 
formation takes the same form as in Eq. (|5.7|) . but with 

G C. The distance d and the unitarity norm are not 
invariant under this transformation, which offers the pos- 
sibility to bring them closer to by a suitable choice of 
the gauge parameters uJa- Hence changes in P[?7, C/^] can 
be quantified by changes in d and the unitarity norm. 

One may wonder why by gauge transformations one 
can improve the result for gauge invariant observables. 
We now answer several aspects of this question. First one 
should remember that the requirement to obtain correct 
expectation values for holomorphic observables does not 
determine the distribution in the complexified space; this 
has been discussed in some detail in the appendix of Ref. 
|25| . In principle there is therefore tremendous freedom to 
select a process having good falloff properties at infinity, 
which can then be sampled efficiently. Gauge invariant 
holomorphic observables allow in addition changing the 
process by using complexified gauge transformations. This 
also changes the process and can be used to prevent it 
from moving out too far along the gauge orbits, which are 
noncompact due to complexification. As will be seen, this 
is in fact necessary in order to keep numerical errors from 
building up and leading to incorrect or unstable results. 

We also note that the procedure of "gauge cooling", 
to be introduced shortly, superficially seems to make the 
process 'less ergodic', because it is suppressing large ex- 
cursions along the gauge orbits in the imaginary direction. 
But here it should be noted that ergodicity in the complex- 
ified configuration space is not necessary for correctness, 
as can be seen already in the simplest one-link U(l) model 
|13) , where the equilibrium configurations contains a delta 
function in the imaginary part. Moreover, ignoring prob- 
lems of numerical precision, the process projected on the 
space of gauge orbits (the quotient space of the full con- 
figuration space by the gauge group) is not affected by the 
cooling, so it is as ergodic as the process without cooling. 

Finally, we note that the formal justification of the CL 
method given in Ref. [13] depends crucially on the holo- 
morphic nature of the drift entering the process. Gauge 
cooling, in the way we implement it, certainly cannot be 
described by such a holomorphic drift. But since we need 
to justify the process only for gauge invariant observables, 
it is only the projection of the process onto the space of 
gauge orbits that matters. The latter, as remarked above, 
is unaffected by gauge cooling if we ignore problems of 
numerical precision. 

As a precursor to gauge cooling, we first discuss a sim- 
ple implementation to illustrate the effect of gauge trans- 
formations on the distribution [T2]. We consider the SU(2) 
model (j5.4p with A'^ = 1 and complex /? = i, and express 
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Fig. 7. SU(2) one-link model: scatter plots of TrU without 
(top) and with (bottom) gauge fixing, a.t j3 — i [12j . 



the matrix U as 

U = cl + iaaba, c^ + baba^l, (5.8) 

where aa are the Pauli matrices. A CL update reads U ^• 
U' = RU, where R was defined in Eq. (|2.20l) . A scatter 
plot of Tr U during the CL evolution is shown in Fig. [7] 
(top). The distribution is very broad and not well- localised: 
indeed the resulting expectation value (TrJ7) disagrees 
with the exact result. However, the model only depends 
on TrC/ = 2c and the full evolution of U contains there- 
fore considerable gauge freedom. This can be constrained 
by diagonalising U after each CL update (gauge fixing). 
The resulting scatter plot. Fig. [7] (bottom), shows a dis- 
tribution which is extremely well-localised, it is effectively 
zero outside the banana shape, and the expectation value 
(Tr U) now agrees with the exact result [12] . We conclude 
that a combination of CL updates and gauge transforma- 
tions has the desired effect. 

The simple gauge fixing used above does not generalise 
to gauge theories in an obvious way. In Ref. |12) it was 
shown that using gauge fixed simulations in the maximal 
axial gauge helps in decreasing the width of the link dis- 
tribution in the imaginary directions. Here, however, we 
use gauge cooling to control the deviation from the uni- 
tary submanifold without specifying a conventional gauge 
condition. This can be achieved by choosing uJa as the gra- 



1 I 




Fig. 8. Gauge cooling in Sh{N, C) brings the links as close 
as possible to S\J{N). The orbit on the left is equivalent to a 
S\J{N) configuration, while the one on the right is not. 

dient of a suitably chosen function, such as the distance d' 
or the unitarity norm. In principle Wq is complex; however, 
choosing Ua (anti)parallel to the gradient yields a unitary 
gauge transformation which is of no interest here and we 
therefore take to be purely imaginary. Starting from 
the distance 

1 1 

<f-]^E]^Tr(c/.f/t-l), (5.9) 

k=l 

we arrive at the following cooling update at site fc, 
Uk ^ Ul ^QkUk, 

Uk-i ^ C/^i - Uk-iOk^ (5-10) 

with 

Qk ^ e-'^'^-f" , (5.11) 

where 

= 2Tr A, (UkUl - Ul^^Uk-i) . (5.12) 

Here e is the Langevin time step and a is a (positive) 
parameter, which can be chosen at will. Below we show 
that choosing a adaptively can lead to very fast evolution 
towards the unitary submanifold. 

In order to demonstrate that a cooling update indeed 
brings the configuration closer to SU(A^), we consider a 
cooling update at site k and find that the distance cf 
changes, to first order in ea, as 

^-d=-^U^)' + 0{{ear). (5.13) 

If the original configuration is gauge-equivalent to a uni- 
tary one, cooling will eventually transform the configura- 
tion into a unitary one. Since — for a unitary config- 
uration, cooling then no longer has an effect. If the con- 
figuration is not gauge-equivalent to SU(iV), cooling will 
bring it as close as possible, i.e. it will minimize d. These 
ideas are sketched in Fig. \E\ 
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Before turning to a numerical solution of CL combined 
with gauge cooling, we first analyse cooling in a one-link 
model analytically. Since we only consider cooling and no 
drift, the form of the action does not come into play; how- 
ever, for definiteness one may consider the models in Eqs. 
(|5.4p and (|5.5p with A'^ = 1. After taking the continuous 
Langevin time limit (e — > 0), the change in rf, to leading 
order in a, is given by 

^=-^fl (5-14) 

with 

fa^2Tr\a{UU^ -U^U) . (5.15) 
Using the Fierz identity, 

A^^ A;;' = 2 [5u5jk - , (5.16) 

this takes the form 

d=-^TYUU^[U,U\ (5.17) 

Specialising to SU(2) to simplify the expression for the 
trace, we find 

cf= -8a(d'2 + 2(l-|c|2)d'+c2 + c*2-2|cp), (5.18) 
where 

c = ^Tr [/, c* = ^Tr U'^ . (5.19) 

To understand this expression, we first note that in the 
one-link model c and c* are invariant under gauge trans- 
formations. Furthermore, for an SU(2) matrix, c = c*. 
Therefore a matrix U e SL(2, C) with a complex trace 
cannot be gauge-equivalent to an SU(2) matrix: d' should 
then remain larger than under cooling. If on the other 
hand the trace is real, the evolution equation simplifies to 

cf=-8a(rf+2-2c2)d', (5.20) 

with the unitary fixed point rf = (the other fixed point 
cannot be reached, since rf -I- 2 > 2c^; here and below 
we exclude the trivial case of the identity matrix, with 
c = c* = 1). From the explicit solution, 

d{t) =c2-l + (l-c2)coth(8a(l-c^)t + const.), (5.21) 

we see that the unitary fixed point is reached exponentially 
fast, 

d{t) - 2 (1 - c^) e-i6a(i-c^)t ^ Q_ (5 22) 

For an SL(2, C) matrix which is not gauge-equivalent to an 
SU(2) matrix, with c* ^ c, the fixed point of the cooling 
evolution is at 

do = |cp - 1 + - c2 - c*2 + |c|4 > 0, (5.23) 

which is again reached exponentially fast. We conclude 
that in the SU(2) one-link model, gauge cooling results 



in an exponential reduction of d — d'o, where d'o is the 
minimum given in Eq. (|5.23p . We note that an exponential 
approach is to be expected on general grounds because 
rf (as well as the other norm) is a smooth function of 
the gauge transformations; hence the gradient will vanish 
linearly at the minimum on the gauge orbit. 

The cooling procedure itself can be improved similarly 
to the dynamics. One method is to use an adaptive step 
size (equivalently, an adaptive cooling parameter a). We 
shall discuss this in connection with the numerical results 
below. Another method is Fourier acceleration [35] ■ To 
demonstrate this, we consider a lattice with sites and 
discuss first the case of U(l) with links 

t/:,,^ = exp[i^:r,M + B^,^,] . (5.24) 

Since the distance as defined above is not bounded by 
for the U(l) case, we use the symmetrised version 

^ = ^ E (t^x,pt/.,M + u-l*u-l - 2) 

a;, /J- 

= ^ (cosh[2B,,^] - 1) > 0. (5.25) 

We consider the gauge transformation 

Bx.^i — > Bx^f_, + y^+f_, - y^. (5.26) 

Let the minimum of d over gauge transformations be at- 
tained for B = and denote the minimal value of d by 
d'o- Then, writing 

Bx,^.^Bl^ + V.,+^-y^, (5.27) 

the gradient of d with respect to the gauge parameters is 

^cf - - 5] ( sinh[2Sj!_^ + 2y,+^ - 2yJ 

- sinh[2B°_^_^ + 2y^ - 2t/,_^]) . (5.28) 

Because the gradient has to vanish if all y^ — 0, B^ satis- 
fies something like a Landau gauge condition: 

E (sinlipi?^,^] - sinh[2i3^^,^]) - (5.29) 

for all X. Linearising Eqs. (j5.28L (I5.29[) both in the gauge 
parameters yx and the fields i?^ ^ we obtain 

-^d = -2 ^ [yx+t. - 2yx + yx-p) = -2Ay, (5.30) 
oyx ^ 

where A is the lattice Laplacean. In this approximation 
thus 

cf=d-o + 2^(grad2/)2 =d'o-2(y,Z\y). (5.31) 

Without Fourier acceleration the gauge cooling equation 
reads, to the same order, 

d^-aiAy,Ay), (5.32) 
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Fig. 9. SU(2) Polyakov chain: evolution of cf without CL up- 
dates, at a fixed gauge cooling parameter a — 0.008, for various 
values of Ni, starting from a gauge transformed unitary chain. 
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Fig. 10. As in Fig. [5J for various values of the gauge cooling 
parameter a and various adaptive implementations, with A'^ = 
1000, on a log-log scale. The dotted line indicates a power 
decay, with power 3/2. 



which gives an exponential approach to the minimum vi^ith 
the rate determined by the softest mode k — 2tt/L. 

To implement Fourier acceleration we calculate the 
Fourier transform of the gradient, divide it by = 
2cosfc^) and transform it back using the inverse Fourier 
transform. This leads to 



d' a{y, Ay) = -a (d - do) , 



(5.33) 



which shows again an exponential approach, but now with 
a rate a. For the nonabelian case we just use this proce- 
dure for every colour, i.e. the variables B^^^^yx acquire 
a colour index a; to leading order in those variables the 
method is justified just as for U(l). 
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Fig. 11. As in Fig. 1101 on a linear-log scale. 



5.1 Gauge cooling for the SU(2) Polyakov chain 

We now turn to the numerical solution of the SU(2) Pol- 
yakov chain (j5.4p with Ng links. We start by considering 
cooling only, i.e. there is no actual CL step in the update. 
Since only the combination ea appears in the cooling up- 
date, we take e = 1. In Fig. [21 we show the dependence 
of the distance rf, see Eq. (15. 9p . on the number of cool- 
ing steps, for a number of chains with Ng links. As initial 
condition we use an SL(2, C) gauge transformed unitary 
chain, i.e. 

;7fe(o) = nuVk^kli^ Vk e su(2), (5.34) 

with fik G SL(2, C) a random gauge transformation. Cool- 
ing is therefore expected to return the chain to being uni- 
tary, with d = 0. This is indeed visible in Fig. [3] and for 
a small number of links this happens exponentially fast 
(note the vertical log scale). The rate, however, decreases 
as the number of links is increased and eventually, for 
Ng = 100, the reduction appears to become nonexponen- 
tial. 

This is further demonstrated in Fig. [TOl where the 
evolution under cooling is shown for a long chain with 
Ng = 1000 links, starting again from a gauge transformed 



unitary chain. The top three lines are obtained as above, 
using three values of a: a larger a results in a faster cool- 
ing, as expected. However, in all cases we note that the 
decay is power-like rather than exponential, with a power 
3/2, as indicated with the dotted line (note the log-log 
scale). It is therefore useful to optimise the reduction of 
d, after a given number of cooling steps, by employing 
an adaptive cooling algorithm. Here the guiding princi- 
ple is to use an effective adaptive aad which is as large 
as possible. However, choosing aad too large may lead to 
instabilities when the cooling drift is large as well, i.e. far 
away from the unitary submanifold. We therefore control 
Q!ad by moderating it appropriately. Explicitly, we have 
implemented adaptive schemes of the form 

a 

"^'^ = mwy ^^-^^^ 

where D{U, depends on the links at the present cool- 
ing step. In Fig. [TU] we present results for D = TrJ7i7^, 
both locally at site k as well as averaged over the chain, 
and D determined by the absolute value of the cooling 
drift, D = Y,ak \fa\INg + Do, with Do a constant. When 
the chain is close to being unitary, both choices reduce to 
a constant effective a, with aad = Oi and a/ Do, respec- 
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Fig. 12. SU(2) Polyakov chain: histogram of the real (top) 
and imaginary (bottom) part of the action, obtained during a 
CL simulation using no cooling (broadest distribution), cooling 
with a constant a (middle) and cooling with an adaptive Qad = 
0.2/(|cooling drift| +0.8) (narrow distribution), using A'^ = 30 
links at /3 = |(1 + Note that the value of a represents 

the product ea. 



tively. From Fig. [TO] we note that the adaptive cooling 
results also contain an interval with sub-exponential de- 
cay, however they are much more effective at the earlier 
stages of coohng. This is demonstrated in Fig. [TTJ where 
the same data is presented, but now using a log-linear 
scale. With an adaptive implementation, after only a few 
cooling steps the distance from SU(2) is greatly reduced, 
drastically changing the distribution in the complexified 
space, as we will see now. 

To proceed, we now include CL evolution with a com- 
plex drift and consider an SU(2) chain with Ni — 30 links 
at /3 = i(l iV3). There is a sign problem due to the 
complex coupling. In order to control the distribution, 
each CL update is followed by a number of gauge cool- 
ing steps. In Fig. [T2]we present distributions for the real 
and imaginary parts of the action S", where the data is col- 
lected after performing a CL update followed by a series 
of gauge cooling updates. Three histograms are shown: 
one without cooling, one with ten cooling steps at a fixed 
a = 0.001 and one with ten cooling steps using an adap- 
tive aad — 0.2/(|cooling driftj + 0.8), where the cooling 
drift is averaged over the chain (note that the value of a 
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Fig. 13. SU(2) Polyakov chain: Langevin time evolution of rf, 
without cooling and using the cooling algorithms of Fig. 1121 
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Fig. 14. SU(2) Polyakov chain: imaginary (top) and real (bot- 
tom) part of the expectation value of the action as a function 
of the number of cooling steps used between CL updates, for 
the two cooling algorithms of Fig. 1121 The dotted lines indicate 
the exact results. 



represents the product ea). We see significant changes in 
the distribution: without cooling the distribution is wide 
and slowly decreasing; with fixed a the distribution is al- 
ready more localised but skirts are still present. Finally 
with the adaptive choice, the distribution is localised very 
well and drops quickly to zero outside its main support. 

These features are also reflected in the distance d from 
the unitary submanifold. This is demonstrated in Fig. [T^ 
where we show the Langevin time evolution of d. Without 
cooling, we find that the average distance equals (d) = 
0.26(2), i.e. relatively far away from the unitary subman- 
ifold. Employing ten steps of cooling with a constant a = 
0.001 reduces this to (rf) = 0.00425(2). Finally with the 
adaptive cooling implementation this is reduced even fur- 
ther to (rf) 0.00029822(4). Cooling therefore brings the 
configuration closer to SU(2). Since the coupling and there- 
fore the action is complex, it should be noted that cf neces- 
sarily has to exceed 0, but apparently it can be very close 
to 0. 

Finally, in order to see whether the CL simulation 
yields a sensible result, we compare the expectation value 
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value as a function of gauge cooling parameter a, using 5 gauge 
cooling steps, for various values of A''^, at fixed /3, k and [17] . 



of the action with the exact result at /3 
obtained by direct numerical integration, 



-0.119121 + ^0.225487. 



(5.36) 



The result is shown in Fig. [T3]as a function of the number 
of cooling steps used between the CL updates. We note 
that the data represent results from simulations using a 
number of Langevin stepsizes and performing an extrapo- 
lation to zero stepsize. For the real part, the result without 
cooling is far off the exact result and not shown. With fixed 
a = 0.001 we observe that the outcome improves when 
the number of cooling steps is increased but there is still 
disagreement after using ten gauge cooling steps. On the 
other hand, using the adaptive implementation, we find 
agreement with the exact results, provided that at least 
six cooling steps are used between CL updates. This agree- 
ment supports our insight based on the requirement of the 
localisation of distributions in the complexified space. 



5.2 Gauge cooling for the SU(3) Polyakov chain 

Since the SU(3) theory at nonzero chemical potential rep- 
resents a genuinely complex problem, we need to verify 
that the understanding gained from the SU(2) case also 
applies to the SU(3) case. We first consider the Polyakov 
chain model with the action (15.51). We refer to 



P 



^Tr ([/,.. 



(5.37) 



as the Polyakov loop. We use again as many as 1024 links 
to simulate the number of variables of a small lattice (the 
results as shown on the smaller chains are, however, al- 
ready representative). Since the procedure is similar to 
the SU(2) case we here concentrate on some results. More 
details can be found in Ref. [TT] . 

In Fig. [15] we compare the expectation value of the 
Polyakov loop (P) with the exact result, for three chains 
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Fig. 16. SU(3) Polyakov chain: histograms of the real part of 
the Polyakov loop, using various combinations of a and the 
number of cooling steps, with Ne — 32 links [171. 
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Fig. 17. SU(3) Polyakov chain: Langevin time evolution of the 
unitarity norm, using various a values {Ne = 32). 



with 16, 24, and 32 links. Between each CL update 5 gauge 
cooling steps are inserted, with different choices for the 
gauge cooling parameter a to enhance the efficiency of 
the cooling. As above, the exact result is reproduced when 
sufficient cooling is applied. As shown in Fig. [THl this con- 
clusion is consistent with results for the Polyakov loop dis- 
tribution for the various a's. Cooling drastically reduces 
the "skirt" of the distribution to a level at which it be- 
comes harmless, as corroborated by the results in Fig. 1151 
We have used both the distance (|5.1D and the unitarity 
norm (|5.2p as the function to be minimised under cooling. 
In Fig. [T7] the Langevin time evolution of the unitarity 
norm, averaged over the chain, is shown. Here we applied 
100 cooling steps between each Langevin update, using 
again various values of a. The cooling dynamics is deter- 
mined by the gradient of the unitarity norm. As above we 
find that substantial cooling stabilises the dynamics and 
keeps the distance from the unitary submanifold under 
control. 

As in the SU(2) case, the evolution of cooling typically 
shows two regimes, a power law behaviour at earlier times 
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Fig. 18. SU(3) Polyakov chain: evolution of the unitary norm 
without CL updates, starting from a gauge transformed uni- 
tary chain, using various cooling implementations {Ne = 32). 



followed by an exponential falloff. We observe the same 
power 3/2 as in the SU(2) Polyakov chain. With a non- 
adaptive a, the power law regime dominates for a long 
time, see Fig. \TE\ (the length of the power law regime in- 
creases rapidly with the number of links). The efficiency 
of cooling can be strongly enhanced by observing that the 
speed with which the minimum is approached is deter- 
mined by the smallest nonzero Fourier mode. Adaptive 
cooling and Fourier acceleration as described above can 
strongly increase the efficiency of cooling, as demonstrated 
in Fig. [THl We conclude therefore that the gauge cooling 
is as efficient as in the SU(2) case and allows one to reach 
the correct results. 



5.3 Gauge cooling for QCD with heavy quarks 

As our final step towards full QCD, we consider QCD at 
nonzero chemical potential with heavy quarks (HQCD). 
This theory retains many of the features of QCD and rep- 
resents an approximation to the latter in the large mass, 
large chemical potential region (heavy dense QCD) [37l 
I38| . For Wilson fermions, with hopping parameter k, it is 
based on a resummed hopping parameter expansion, re- 
taining only time-like Polyakov loops, which can formally 
be obtained from the limit k — >■ 0, /i — >■ oo, with kc'^ 
fixed. This approximation, and variants thereof, have re- 
peatedly been used to obtain information about the phase 
structure of QCD in the /x — T plane [55115^ . Here we use 
a "symmetrized" version [5] of the model, with the action 

S' = 5YM-lndetM(^). (5.38) 

Here Sym is the standard Wilson action on a x A^^ 
lattice with gauge coupling j3, and the determinant takes 
the form 

det M = [| det (l + hd'/^V^ ^ det (l + /le^^/^-p^ ^ , 

(5.39) 



100 




5 10 15 20 25 30 



Langevin time 

Fig. 19. QCD with heavy quarks (HQCD): evolution of the 
distance d with and without gauge cooling, at /i = and 1.3, 
for /? = 5.9 on a 8^ X 6 lattice, k = 0.12 throughout [17] . 

where T = I/At, lattice units are used throughout, h = 
(2k)^^, and the (conjugate) Polyakov loops are 

- n '^(-.x),4, V^^V^ = 1. (5.40) 

T=0 

Below we also consider the traced Polyakov loops, P^'^^ — 
■itrT'x The factor containing e~^l'^ is of course neg- 
ligible in the heavy dense limit, but it is required for the 
symmetry [det M(/i)]* = det M(-/i*). 

A detailed analysis can be found in Ref. JLZJ: here we 
discuss a selection of results, some of which are new. As 
shown in Fig. [191 the effect of gauge cooling is as important 
as in the Polyakov chain model. In this figure we present 
the Langevin time evolution of the (symmetrised version 
of the) distance [T7], for a number of different cases. When 
a = 0, no cooling is applied. The unitary submanifold 
is unstable and the distance increases exponentially. At 
nonzero /i this happens immediately, while at vanishing [i 
it takes some time for the instability to set in. At vanish- 
ing /z this growth can of course easily be avoided by an 
occasional re-unitarisation of the links; however, this pos- 
sibility is not present at nonzero /i and therefore does not 
provide a valid solution. With cooling (a > 0; we used 12 
gauge cooling steps between each CL update) , we observe 
that the distance remains under control and in the case of 
/i = 0, the dynamics remains on the unitary submanifold 
(within numerical precision). 

The baryon density, normalized with the density at 
saturation (risat = 6 in this model), is shown in Fig. [20] 
as a function of [i. We observe an onset, followed by a 
rapid rise to saturation. We note that the entire region 
from /i = to far into the saturation regime is covered. 
Also shown is the average phase factor in the full theory, 
defined by 

, / detM(;.) \ 
\ I \detM(-^)/- ^ ' 

This observable shows that between onset and saturation 
the sign problem is severe, such that other methods such 
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Fig. 20. HQCD: baryon density and average phase factor as 
a function of ^ for /? = 5.9 on a 8^ X 6 lattice ITfl. 
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Fig. 21. HQCD: spatial and temporal plaquettes (top) and 
Polyakov loops (P) and {P~^) (bottom) as a function of fi: 
comparison between CL and reweighting at /3 = 5.9 on a 6^ 
lattice, using a = 1 with 12 gauge cooling steps. Large errors 
at large fi affect RW only [TT] , 



as reweighting will be ineffective here. Nevertheless, CL 
appears to work well. We note that the average phase 
factor is a gauge invariant observable, just as the density 
or the plaquettes, presented below. The role of cooling in 
improving results for gauge invariant observables has been 
discussed at some length below Eq. (|5.7p . We note that 
without (sufficient) cooling, observables do not converge 
to the expected values, similar as in Figs. [U and [12] for 
the Polyakov chain models. 

For QCD with heavy quarks, no exact results or re- 
sults obtained with sign-problem free approaches are avail- 
able. We therefore compare expectation values of (traced) 
Polyakov loops and temporal and spatial plaquettes with a 
refined reweighting (RW) analysis as described in Ref. . 
In Fig.[5T]we show the spatial and temporal plaquettes as 
well as the (conjugate) Polyakov loops. Very good agree- 
ment is found in the region where reweighting is feasible; 
for ^ > 1 the RW data rapidly deteriorate, due to the sign 
problem. 

In Fig. [22] we demonstrate that it is possible with CL 
to go from the deconfined to the confined phase by varying 
/3 at fixed /j,. We observe that the Polyakov loop averages 
agree with the RW results for all /3 values considered. For 
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Fig. 22. HQCD: Polyakov loops (P) and (P"^) (top) and spa- 
tial plaquettes (bottom) as a function of /?: comparison between 
CL and reweighting at /.i = 0.85 on a 6* lattice [T7]. 
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Fig. 23. HQCD: temporal and spatial plaquettes as a function 
of /3: comparison between CL and reweighting at ^ = 0.85 on 
a 10^ lattice. 



the plaquette expectation values we note that they agree 
except for smaller /3 values. In Fig. [53] we show results on 
a larger 10''' lattice: increasing the lattice size does not ap- 
pear to move the threshold in (3 below which the plaquette 
results deteriorate. First indications are that on the 10* 
lattice the transition to the confinement region takes place 
at /3 ~ 6. This suggests that if a deterioration of the CL 
results occurs, it is a small (fixed) (3 effect and it is possi- 
ble to go deep into the confining region by increasing the 
lattice size, a situation very different from the one in the 
XY model. If this conclusion remains true in full QCD it 
means that we can safely approach the continuum limit of 
the latter in both the confinement and the deconfinement 
phases, which is an exciting prospect. 



6 Summary 

Complex Langevin dynamics can solve the numerical sign 
problem appearing in lattice field theories with a com- 
plex action or Boltzmann weight. However, success is not 
guaranteed. In this paper we reviewed the reasons why a 
simple proof, valid in the case of real Langevin dynamics, 
is not available and outlined an alternative route. We indi- 
cated how problems can manifest themselves in numerical 
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simulations. Particular emphasis is put on the properties 
of the positive probability distribution in the complexified 
space, which is effectively sampled during a CL process. 
We argued that a proper understanding of this distribu- 
tion can be used to provide support for the validity of CL 
results - or discard the results as a failure. These insights 
are demonstrated by contrasting results obtained in two 
three-dimensional models with a nonzero chemical poten- 
tial: the XY model and the SU(3) spin model. 

In the second part we applied the lessons learnt to 
theories with SU(iV) gauge symmetry. It is demonstrated 
how the freedom to perform SL(iV, C) gauge transforma- 
tions can be used in a constructive manner to control the 
distribution in the extended space. By interspersing CL 
updates with SL(A^, C) gauge cooling steps, the distance 
from the unitary submanifold can be reduced. We have 
demonstrated that this can cure the problem of conver- 
gence to the wrong result. Adaptive cooling proves to be 
useful in that it can lead to a surprisingly quick reduction 
of the deviation from the unitary submanifold, thereby 
greatly improving the convergence properties of CL. 

The gauge cooling techniques discussed here build on 
the results presented in Ref. |17j . in which they were ap- 
plied to QCD with heavy quarks. We are currently ex- 
tending the analysis, taking into consideration the insights 
gained here, and hope to report on our findings in QCD 
in the near future. 

We thank Simon Hands for the opportunity to present this 
work. GA and LB are supported by STFC. ES and lOS were 
partially supported by Deutsche Forschungsgemeinschaft. 
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